Pattern formation in self-propelled particles with density-dependent motility 
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We study the behaviour of interacting self-propelled particles, whose self-propulsion speed decreases 
with their local density. By combining direct simulations of the microscopic model with an analysis of the 
hydrodynamic equations obtained by explicitly coarse graining the model, we show that interactions lead 
generically to the formation of a host of patterns, including moving clumps, active lanes and asters. This 
general mechanism could explain many of the patterns seen in recent experiments and simulations. 
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Collections of self-propelled (SP) particles provide the 
most common realization of active matter, the study of 
which constitutes a rapidly growing area of research in 
physics [1 1. Examples of SP particles are bacteria, cells J2) 
and actin filaments "walking" on a carpet of immobilized 
molecular motors 151 . 

The term "active" is used to contrast these systems with 
their passive counterparts, such as solutions of diffusing 
Brownian particles. Active systems exhibit a much richer 
physics than their passive counterparts. Most important 
for us, they have a far larger tendency to form patterns. 
For instance, bacterial colonies of e.g. E. coli or S. ty- 
phimurium growing in the lab can self-organize into crys- 
talline or amorphous arrangements of high-density bacte- 
rial clumps [|4), while biofilms form even more elaborate 
patterns such as microbial honeycombs, essentially hexag- 
onal lattices of low-density spots, or voids [5|. Similarly, 
actin in high density motility assays [ 3 1 organize in moving 
spots or stripes as well as traveling waves. 

What is the mechanism underlying the formation of 
these "active patterns"? One may expect that, as the un- 
derlying constituents of each system are so different, the 
answer to this question should also be system-specific. If 
we are to capture all details of a given active pattern, this 
is indeed likely to be the case. Yet, a fascinating possibil- 
ity is that there may exist some generic origin of many of 
these patterns, stemming from a few universal key features 
of activity, linked to its inherent non-equilibrium nature. 
In some cases, pursuing such minimal descriptions can be 
very rewarding. A well-known example is the hydrody- 
namic theory of flocking proposed by Toner and Tu in JSJ, 
which was inspired by the "agent-based" model of Vicsek 
et al. 1UG1. The latter studied the dynamics of an ensemble 
of SP particles subjected to aligning interactions, whose ul- 
timate origin may be hydrodynamic or collision-dominated 
in the cases of bacteria and actin filaments, or more com- 
plex for bird flocks or fish schools. Universal features suc- 
cessfully predicted by generic flocking models are spon- 
taneous motion, giant density fluctuations, and the emer- 
gence of complex spatiotemporal active patterns HOIS. 

The original Vicsek model considers point particles of 



fixed speed and includes no interactions between the SP 
particles other than a rule that aligns their velocities. Re- 
cently, focus has shifted onto specific models where addi- 
tional interactions are included, most commonly steric re- 
pulsion fTOttToll . Our aim here is to develop a more generic 
model for interacting SP particles. Interactions are incor- 
porated in our model by assuming that the motility of the 
SP particles is a decreasing function of their local den- 
sity ifTTl . One may envisage several physical mechanisms 
responsible for a decay of the propulsion velocity with den- 
sity: here we highlight just two. First, such a slowdown 
may arise due to local crowding and steric hindrance, just 
as in IfTTl [T2l [141 [T31 . An alternative mechanism can be 
provided by biochemical signaling such as quorum sensing 
in bacterial colonies, as recently explored theoretically lfT8l 
and experimentally [19]. This second mechanism may lead 
to slowdown even in dilute suspensions. Our work de- 
scribes the results of simulations of a microscopic SP par- 
ticles model with both interactions and alignment rule, the 
derivation of the corresponding hydrodynamic description 
of the model in terms of a density and a polarization field, 
and an analysis of the continuum theory. It therefore pro- 
vides a direct bridge between microscopic and continuum 
models, which allows us to identify universal mechanisms 
driving pattern formation in interacting SP particles. As we 
shall see, interactions lead to an even larger repertoire of 
patterns in active particle suspensions than obtained in con- 
ventional Vicsek models. These include moving clumps, 
lanes and asters, and qualitatively match the patterns found 
experimentally, e.g. in J3J. Our hydrodynamic equations 
also provide an understanding of the origin of the various 
patterns and allow us a clear comparison with other inter- 
acting active particle models. 



We consider a modified version of the Vicsek model Q, 
where the particles interact via a pairwise aligning forcing, 
which simplifies the coarse graining of the microscopic 
model. In 2D the position r.; and direction, identified by 
an angle 9i (or a vector e e .), of the ith particle evolve ac- 
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cording to 

(1) 

where 7 and e are parameters describing the strength of 
alignment and fluctuations respectively, and fj(t) is a Gaus- 
sian white noise with zero mean and unit variance. F con- 
trols the alignment interactions between the spins. For sim- 
plicity, we choose F(9, r) = sm(6)/irR 2 if |r| < R and 
otherwise, though its precise shape does not dramatically 
affect the physics. In the v — > limit, our model is an 
off-lattice analogue of the XY model for a ferromagnet, 
hence we call it the flying XY model. This differs from 
other models of SP particles where alignment is explicitly 
due to 'collisions', and in which the interaction strength 
vanishes as v — > ll20l . Such cases can however be recov- 
ered by taking 7 oc v. Last, a density-dependent velocity 
is introduced in the model by stipulating that v depends 
on the number of particles n within a given radius R n , as 
v(n) = v Q e~ Xn + Vi, where v 3> V\ > are the dilute 
and crowded limiting velocities respectively, and A > 
controls the decay of the motility decreases with increas- 
ing density. Hereupon we restrict to R n = R. 

Fig. shows a representative phase diagram of the fly- 
ing XY model in the (e, A) plane. For small A, when v is 
practically constant, the phases observed are the same as 
those in the literature on flocking models HUT)- Namely, 
at high e we find a disordered, homogeneous state (region c 
in Fig. 1 A), followed by a polarly ordered phase with high 
density stripes (named stripy phase, b in Fig. 1A) below a 
critical noise value. For even lower e, we observe a 'fluc- 
tuating flocking state' (region a) with polar order and large 
density fluctuations - this state is very close to the one de- 
scribed in Ref. [9| and we do not discuss it further here. 
These phases are expected - as we shall see below the hy- 
drodynamic equations for our model map to those for the 
Vicsek model J6j|9) when v is constant. 

Above a critical value A c (e), new patterns appear. Due to 
the density-dependent motility, the SP particles cluster via 
a self-trapping mechanism through which they assemble 
and slow-down, creating a positive feedback loop akin to 
the one discussed for models of bacteria |[T8ll . This process 
leads to the formation of high density clumps which slowly 
coarsen towards a fully phase separated steady state. The 
Vicsek-like alignment tendency greatly affects this instabil- 
ity. On one hand, the critical value A c (e) decreases almost 
to zero with decreasing e. Furthermore, the presence of po- 
lar order promoted by the alignment changes the nature of 
the clusters. In Fig. 1A we identify at least three distinct 
patterns, of which snapshots are shown in Fig. IB. When e 
is small, rather than structureless dots, the clusters show an 
orientational order and move coherently: they form "mov- 
ing clumps" (pattern d(i) in Fig. 1). For low e and large 
A the moving clumps merge into bands, or lanes (labelled 
as d(ii)) - within these, however, particles move parallel 
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FIG. 1: (color online) (A) Phase diagram in the (e, A) plane, for 
N = 3000, L = 10, 7 = 0.16, v = 2 and vi = 0.1. Blue filled 
circles on the phase boundary correspond to peaks in the variance 
of the particle density, while green squares separate states with 
zero mean orientation from states with nonzero mean orientation. 
Phases are labelled as per discussion in the text. Horizontal and 
vertical red lines indicate linear instabilities towards clustering 
and ordering, respectively. (B) Snapshots of the stripy (b), aster 
(e), moving clumps (d(i)), lane (d(ii)) patterns. The crosses in 
A correspond to the snapshots in B. Particles are color coded by 
direction, with blue horizontal and red vertical. 

rather than perpendicular to the band, in contrast with the 
A — > stripy phase. This is reminiscent of the 'streaks' 
of actin filaments observed in [3|. In the disordered, high 
e phase, the clusters instead diffuse randomly, and are on 
average stationary. 

Here, a temporal average of the particle orientation pat- 
terns shows that the clusters are asters (the aster phase is 
labelled as e in Fig. 1). However, as discussed in greater 
detail below, the orientation in the aster is non-standard: 
particles point towards its center at the core, but they co- 
herently point outwards in its periphery. We stress that 
moving clusters, lanes and asters are not observed either 
in the standard Vicsek model, or in its standard mean field 
continuum description used in Ref. [9|. They are a conse- 
quence of the interplay between self-propulsion, alignment 
and density-dependent motility: switching off any of these 
ingredients would thus greatly reduce the pattern forming 
potential of the model. 

To get a better understanding of the pattern formation 
process, we now discuss how to coarse grain the micro- 
scopic dynamics ([TJ to obtain a macroscopic description of 
the model. They are two candidates for the hydrodynamic 
fields: the particle density p, which is conserved, and the 
local alignment, or polarization, vector P, on symmetry 
ground. Note that "hydrodynamic" here means slowly 
varying in space and time - the dynamics of the underlying 
fluid is not included in our modeling. Following Ref. J22), 
we start with the microscopic Eq. ([T]l and use I to calcu- 
lus |[T8l to write down a stochastic dynamical equation for 
the evolution of /(r, 9, t), the density of particles at posi- 
tion r with angle 8, as follows 

d t f(r, e)+e e . V[vf] = e^l - ^V^fv 
-!§- e f de'dr'f(r',9')f(r,e)F(O'-0,r-r'). 

(2) 



3 



The second term on the left hand side describes familiar 
advection, but with one important difference: the velocity 
v appears inside the gradient. This is precisely what leads 
to the instabilities responsible for the new patterns in the 
simulations. By dropping the conserved noise ^2efrj, we 
obtain a non-fluctuating kinetic equations for the flying XY 
model. Following Berlin et al. l20ll . we Fourier transforms 

Eq.t2lto get equations of motion for f k = J fe lk6 d8. Us- 

ike 



ing^r/(#) = Z k f k e- lke and 2ttF{0) 
we obtain a hierarchy of equations: 
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where all sums run from — oo to +oo. In principle, F is 
slightly non-local in space so that the second term of the 
r.h.s. of Eq. Q should retain a spatial integral. We are 
however interested in the hydrodynamic, large-scale, de- 
scription of the system, a limit in which R is very small 
and we assume F to be perfectly local ll23l . To obtain 
mean field equations for the hydrodynamic variables, we 
note that f k for k = is simply the density, p, whereas the 
real and imaginary part of fi are the x and y component of 
pP respectively (as we work in 2D we can identify vectors 
with complex numbers). By writing out in full the k = 
case of Eq. Q, we then find that the density field obeys the 
continuity equation 

d t p = -V • (vW), (4) 

where W = pP. To make further progress, we now 
assume that we are not too deeply in the ordered phase, 
so that f(9) is to first order approximation homogeneous, 
hence higher Fourier components (f k for k > 3) may be 
neglected. Following EOl . we further assume that f 2 is a 
fast variable, so that f 2 — 0. After lengthy but straightfor- 
ward algebra, we obtain the following equation for W, 
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-tW(V • W) - — v(W ■ V)W + C(V 2 ) 

The second term on the l.h.s. of Eq. Q describes self- 
advection of particles and breaks Galilean invariance 0. 
The first two terms on the right-hand side describe the stan- 
dard spontaneous symmetry breaking leading to polar or- 
der and flocking for sufficiently small e in the Vicsek model 
at A = 0. The third, pressure-like term, — |V(wp), is 
the most relevant one in our work, as it is responsible for 
the clustering instability observed in Fig. 1 when A ^ 0. 
Higher order terms in V and W have minor effects on pat- 
terns and will be discussed elsewhere. 
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FIG. 2: (color online) (A) Phase boundary for the flying XY 
model when A = 0, showing the critical value of e as a func- 
tion of 7. Blue points for v = 2.0, red for v — 0.5. Inset: data for 
v = 2.0 for smaller values of 7, showing good agreement. (B) 
Phase boundary for e = 5, 7 = 0.16. In all cases L = 10 and 
N = 1000. 



Having written down the mean field equations of mo- 
tion, Eqs. Q and ([5]), we can now assess how their predic- 
tions compare with the full simulations of the microscopic 
model. The continuum theory predicts an order-disorder 
transition at e c = |7Po- F° r e > e c there is a stable ho- 
mogeneous disordered state, with p = p and W = 0. 
For e < e c the equations yield a homogeneous ordered or 
flocking state with p = p Q and W = W Q i, where we 
have chosen the x axis along the direction of broken sym- 
metry and Wq = A/8e(e c — e)/7 2 - The mean-field tran- 
sition at e c does not depend on A and coincides with that 
of the equilibrium XY model. The order-disorder phase 
boundary predicted by the theory is compared to its nu- 
merical counterpart in Fig. 2A. We then study the linear 
stability of the homogeneous disordered state at e > e c 
against spatially inhomogeneous fluctuations. It is straight- 
forward to show that when A 7^ the homogeneous disor- 
dered phase becomes unstable for all wavenumbers when 
v (Po) + Po v '(Po) < 0. This instability, referred to as a 
clustering instability, arises due to the term - |\7(fp) in 
the equation for W. The threshold between homogeneous 
and clustered phases found numerically at large e is close 
to but below the prediction (Fig. 2B). This is reasonable, as 
the linear stability can only access the spinodal line: fluc- 
tuations may trigger phase separation for lower A. 

To go beyond the simple linear stability analysis of the 
homogeneous disordered state, account for the effect of the 
non-linear terms, and hence explore the range of patterns 
compatible with our hydrodynamics equations, we solved 
Eqs. Q and Q numerically, by means of a standard fi- 
nite difference scheme. In order to enhance the stability 
of our algorithm, we included a diffusive term DX7 2 p on 
the right hand side of Eq. Q. Our numerical results show 
that all the five patterns, or phases, observed in the mi- 
croscopic simulations (fluctuating flocking state, moving 
stripes and lanes, static asters and moving clumps) can be 
found within Eqs. Q and ^ - Fig. 3 portrays a com- 
parison of the A / patterns. Interestingly the origin of 
the atypical asters can be directly read in Eq. Q. In the 
steady-state, low gradient, small W approximation, Eq. Q 
reduces to (7/9/2 — e)W = X7(pv/2) and V(vp) thus acts 
as an ordering field for W. Along the radius of an aster, 
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FIG. 3: (color online) Patterns found for A / in the microscopic 
simulations (left column) and in the numerical solution of the 
hydrodynamic equations (right). Tables show dimensionless pa- 
rameter values: A = XpQ, v~i = vi/vq, 7 = -ypo/e, D = De/vq. 
Arrows show the W field, colors the density (red: high; blue: 
low). 



the density increases towards the center whereas the veloc- 
ity decreases. Their product can thus be non-monotonous, 
which makes W change direction, whence the atypical 
asters seen in the microscopic simulations. In the contin- 
uum simulations, even though \7pv can change sign, the 
presence of the diffusive terms disallows sharp gradients 
and we did not find parameters for which \/pv was domi- 
nating. We could, however, end up with both inward point- 
ing or outward pointing asters, corresponding to phases 
with high-density clumps (at small A, shown in Fig. 3C) 
or low-density voids (at larger A, similar to those discussed 
in 01, not shown). 

We have shown that a density-dependent motility in our 
flying XY model, which is a close relative of the Vicsek 
model, yields new patterns in suspensions of SP particles. 
Such patterns include moving clumps, lanes, and asters, as- 
sociated with high-density clumps, or voids. All these pat- 
terns have experimental counterparts [ 3-5 1 . By explicitly 
linking the microscopic and coarse grained mean field dy- 
namics, we were able to identify the key ingredients which 
trigger the appearance of the new patterns in the "pressure 
term" — |V(u / o): when this turns negative, new patterns 
may form. Importantly, the patterns we see are not very 
sensitive to the precise form of v(p). For instance, steric 
hindrance results in velocities that typically decrease lin- 
early with density 11241 and would give similar instabilities. 

We close with a comparison with other models fea- 



turing patterns similar to ours. Continuum models for 
microtubule-kinesin solutions leading to aster formation 
have been proposed in [25 1. The resulting equations of 
motion for the microtubule polarisation included a phe- 
nomenological term of the form 5V(p) with S > 0, and p 
the density of motors bound to microtubules. This is qual- 
itatively similar to our equations with negative — |V(up). 
In the A = limit, Refs. J9j|20l show that asters are absent 
if the prefactors in the non-linear terms in the continuum 
equations are obtained via a systematic coarse-graining of 
a system of Vicsek SP particles or SP hard rods. Ref. 11261 
shows, however, that asters do appear if these prefactors 
are tuned as independent parameters, although in this case 
the asters have fixed polarity. 

Finally, Peruani et al. lfT5l studied a microscopic lattice 
variant of the Vicsek model, and also found asters and mov- 
ing clumps, dubbed traffic jams and gliders respectively. 
The underlying passive model in Ref. |[T5l is the 4-state 
Potts model, which is in a different universality class than 
the XY model, to which our dynamics reduces in the pas- 
sive limit. However, the active patterns are similar in the 
two models. This is again naturally explained by our the- 
ory, as their origin in |[T5l lies in the slowdown of parti- 
cles due to crowding jamming, which brings up an effective 
"pressure term" analogous to the one in Eq. Q. A density- 
dependent motility, induced either by steric hindrance or by 
crosslinkers between actin filaments, may also at the basis 
of the formation of similar patterns in the actin-walker ex- 
periments in J3). 
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